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Abstract. We describe a powerful methodology for numerical solution of 
3-D self-gravitational hydrodynamics problems with extremely high resolu- 
tion. Our method utilizes the technique of local adaptive mesh refinement 
(AMR), employing multiple grids at multiple levels of resolution. These 
grids are automatically and dynamically added and removed as necessary 
to maintain adequate resolution. This technology allows for the solution of 
problems in a manner that is both more efficient and more versatile than 
other fixed and variable resolution methods. The application of AMR to 
simulate the collapse and fragmentation of a molecular cloud, a key step 
in star formation, is discussed. Such simulations involve many orders of 
magnitude of variation in length scale as fragments form. In this paper we 
briefly describe the methodology and present an illustrative application for 
nonisothermal cloud collapse. We describe the numerical Jeans condition, 
a criterion for stability of self- gravitational hydrodynamics problems. We 
show the first well-resolved nonisothermal evolutionary sequence beginning 
with a perturbed dense molecular cloud core that leads to the formation 
of a binary system consisting of protostellar cores surrounded by distinct 
protostellar disks. The scale of the disks, of order 100 AU, is consistent with 
observations of gaseous disks surrounding single T-Tauri stars and debris 
disks surrounding systems such as /? Pictoris. 
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1. INTRODUCTION 

Star formation ultimately involves the collapse of an interstellar molecular 
cloud core from an initial density oi p < 10~^^ g cm~^ and size ~ 10^^ cm to 
a final young stellar object of density p > lO"-*^ g cm~^ and size ~ 10^^ cm. 
The collapse of a portion of a cloud may lead to fragmentation, which is 
crucial to establishing key parameters of binary stars: the distributions of 
mass ratios, periods, and orbital eccentricities. It may also be essential to 
the formation of larger groups of stars and to the determination of the stel- 
lar initial mass function. However, this enormous dynamic range presents 
a formidable obstacle to obtaining an accurate numerical solution, as the 
flow must remain well-resolved throughout the evolution. Fixed-resolution 
methods cannot be used to simulate such a 3-D collapse in a practical 
amount of time using current computers. 

In order to treat the enormous dynamic range presented by gravita- 
tional collapse calculations and address the problems associated with alter- 
native methods such as smoothed particle hydrodynamics (SPH) and static 
nested-grid codes (Burkert & Bodenehimer, 1993), we have developed a 3-D 
self-gravitational adaptive mesh refinement (AMR) code. Specifically, SPH 
suffers from two major drawbacks. First, it is typically very diffusive, and 
has great difficulty tracking shocks and contact discontinuities accurately 
in multi-dimensional flows. Second, although it is an adaptive method, it 
is only adaptive in the sense that the local smoothing length h tracks mass 
in a Lagrangian sense, i.e. in 3-D, h oc p~^/^. Furthermore, in the context 
of isothermal collapse, in which the Jeans length Aj scales as Aj oc p~^^'^^ 
the effective resolution, specified by the ratio \j/h varies as p~^/^ - the ef- 
fective resolution of the calculation is degraded as collapse proceeds, and if 
the collapse is not arrested first, the calculation will eventually become un- 
derresolved. Given that in SPH one cannot choose a priori which regions of 
the flow need further refinement, one's only recourse to increasing the res- 
olution of a calculation is to increase the total number of particles. Lastly, 
static nested-grid codes suffer from the drawback of imposing a static set of 
grids from the beginning of the calculation, and hence are unable to resolve 
arbitrarily-located fragments. 

Our AMR method dynamically resizes and repositions grids and inserts 
new, finer ones within them according to adjustable refinement criteria 
(unlike the static nested grid scheme of Burkert & Bodenehimer 1993). 
We have applied our method to great effect; prior to the addition of self- 
gravity, we have used AMR to study astrophysical problems including the 
interaction of supernova blastwaves with interstellar clouds (Klein, McKee, 
k Colella 1994; Klein, & McKee 1994; Klein, McKee, k Woods 1995), X-ray 
heated coronae and winds from accretion disks (Woods et al., 1996), and the 
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collision of interstellar clouds (Klein k. Woods, 1998) with unprecedented 
high resolution. 

We present the methodology behind our 3-D AMR self-gravity code 
and describe application of our work to the collapse and fragmentation of 
molecular clouds. Our discussion closely follows Truelove et al. 1998 and 
Klein 1998. 

2. METHODOLOGY 

2.1. SELF-GRAVITATIONAL HYDRODYNAMICS 

The basic governing equations of our model arc the Eiilcr equations of 
hydrodynamics in 3-D, including effects of self-gravitation, and Poisson's 
equation for the gravitational potential. We term these the gravitational hy- 
drodynamics (GHD) equations. The total non-gravitational energy per unit 
mass, E, is related to the internal energy per unit mass, e, hy E = e + ^v^. 
In turn, P is defined by an equation of state. In our isothermal work, we 
adopt the ideal gas law P = (7 — l)pe. In recent work (Fisher, Klein, & Mc- 
Kee 1998) we have used a nonisothermal EOS that combines isothermal and 
nonisothermal components such that P{p) = (?gP + Kp^ where 7 = 5/3 and 
K is chosen so that the isothermal and adiabatic components balance at a 
critical density determined by detailed 1-D radiation-hydrodynamic calcu- 
lations of Masunaga, Miyama, & Inutsuka (1998). A typical critical density 
for this transition is 5 x 10~^"^ g cm~^. As written above, there are eight 
equations in eight unknowns. We assume periodic boundary conditions on 
these equations. 

Our code methodology consists of three components to efficiently solve 
this system. The first component is a hyperbolic solver that employs an 
implementation of the Godunov method (see Colella & Woodward 1984 
and LeVeque 1992) for solution of the Euler equations of gas dynamics. 
The second major component of our code methodology is an elliptic solver 
that utilizes an AMR multigrid method to solve Poisson's equation. These 
elements operate within the third component, an adaptive mesh refinement 
framework. 

2.2. HYPERBOLIC SOLVER 

2.2.1. Overview of Godunov Implem,entation 

To solve the full 3-D system of hydrodynamics equations, we use an operator- 
split method in which we solve in each coordinate direction independently 
in a cyclic sequence. In 3-D, this operator splitting can be written schemat- 
ically as 

^ L,L,L,[L,L,L,(Q")], (1) 
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where the superscript n indicates the n time level, Q" indicates the state 
vector at the n*^ time level, and Li is the update operator in the i^^ coor- 
dinate direction based on the higher-order Godunov method described by 
Colella & Woodward (1984) and Bell, Colella, k Trangenstein (1989). That 
is, each three-dimensional update step from to t^'^^ consists of 3 one- 
dimensional update steps, one per coordinate direction. At the end of the 
cycle from t" to the solution is second-order accurate. The interested 

reader can refer to Truelove et al. (1998) and Klein (1998) for the details 
of the scheme. 

2.2.2. AMR Considerations and AMR Multigrid 

When solving a hyperbolic system on a hierarchy of computational grids, 
care must be taken that the solutions on finer grids are reflected on the 
neighboring and underlying coarser grids. Our code uses refluxing and av- 
eraging down procedures described in detail by Berger & Colella (1989) 
to update the solutions on these coarser grids to account for solutions on 
the finer grids. In refluxing, the solution in a coarse cell adjacent to a fine 
grid is updated after the fine grids have been advanced in time by adding 
the effects of a differential flux acting over At at its face adjacent to the 
refined region. This differential flux is equal to the difference between the 
flux at this face as computed on the coarse grid and the sum of the fluxes 
at this face as computed on the fine grid. The solution in a refined coarse 
cell is simply overwritten by averaging down the solutions on the fine cells 
it contains. 

In contrast to the hyperbolic solution, in which only neighboring grids 
communicate at each time step, the solution of Poisson's equation for the 
gravitational potential is a much more tightly coupled process. The multi- 
grid method is a natural choice to use as an element of the elliptic solver 
in an AMR calculation (Truelove et al. 1998). 

2.3. AMR FRAMEWORK AND THE REFINEMENT CRITERION 

As discussed above, the basic hydrodynamic method we use is a higher- 
order extension of Godunov's method of a type discussed in Colella &: 
Woodward (1984). This algorithm is second-order accurate for smooth flow 
problems, and has a robust and accurate treatment of discontinuities (Klein, 
McKee, & Colella, 1994). We further enhance the efficiency and accuracy 
of our calculations by using local adaptive mesh refinement. 

The version of the algorithm used for the present work is rather elabo- 
rate; a detailed description is given in Berger &; Colella 1989. 

AMR contains five separate components. The error estimator is used 
to estimate local truncation error and will be described subsequently. The 
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grid generator creates fine grid patches that cover the regions that need 
refinement. Data structure routines manage the grid hierarchy allowing ac- 
cess to the individual patches. Interpolation routines initialize a solution 
on a newly created fine grid and also provide the boundary conditions for 
integrating the fine grids. Flux correction routines ensure conservation at 
grid interfaces by modifying the coarse grid solution for coarse cells that 
are adjacent to a fine grid. 

In AMR, the computational volume consists of a hierarchical grid struc- 
ture. A base Level grid fills the computational volume, discretizing it on a 
rectangular grid with a resolution of Axq in each direction. Multiple Level 
1 grids of finer resolution A.^! = Axo/ri may be embedded within it, where 
ri = 4 is a typical choice. Higher levels may be defined embedded in under- 
lying levels in a similar fashion. Grids at Level L always span an integral 
number of cells at Level L — 1, i.e., partial cell refinement is not permitted. 
FTirthermore, an interior grid at Level L is always nested within a grid at 
Level L — 1 such that there is a buffer region of Level L — 1 cells surround- 
ing it. In other words, a grid at Level L within a grid at Level L — 1 never 
shares a boundary with the Level L — 1 grid, unless it shares a boundary 
with the computational domain. 

An integration step on this adaptive grid proceeds as follows. We are 
using a single timestep for all grids, determined from the most stringent 
timestep over the entire physical domain. The updating of the data on 
the locally refined grid structure is organized around the grouping of cells 
into rectangular grid patches, each typically containing hundred to several 
thousand grid cells. The AMR code passes a grid patch to an integration 
routine which returns the updated cell values as well as the edge-centered 
conservative fiuxes used in the update. The overheads in both CPU and 
memory associated with the adaptive mesh structure have been kept quite 
small, relative to irregular grid schemes. Typically, 80% - 90% of the total 
execution time is spent advancing cells in time using the hydrodynamic and 
elliptic solvers, while the memory required is that needed to store two copies 
of the solution on all of the grids. These overheads are low because they 
are determined by the number of rectangles into which the AMR solution 
has been divided; as opposed to being determined by the number of grid 
cells, as is the case with irregular grid adaptive algorithms. 

A key component of an AMR code is the grid generator. In our code, 
this procedure is broken into two steps. In the first step, a specified property 
is measured in each cell, and cells requiring refinement are flagged. In the 
second step, the distribution of flagged cells is analyzed to determine the 
number, sizes, and locations of grids to be inserted at the next finer level of 
resolution. These finer grids will always include every cell that was flagged 
for refinement, but they may also include additional cells that were not 
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flagged. The degree to which the reflnement is concentrated in the ceUs 
that require it is termed the grid efficiency. The grid efficiency is minimal 
when the smallest rectangular solid containing all flagged cells is refined. 
The grid efficiency is maximal when the only cells refined are those that 
were flagged. 

Our general process for deciding upon a suitable level of reflnement de- 
pends on a method for estimating the local truncation error of the Godunov 
scheme; this error estimator determines where the solution accuracy is in- 
sufficient. We estimate truncation error using Richardson extrapolation; we 
take data on the grid where the error is being estimated and evolve two 
timesteps. This solution is then combined with the data on a grid spatially 
coarsened by a factor of two and integrated for one timestep. The local 
truncation error can be shown to be proportional to the normed functional 
that is the difference of these two operations. 

In the nonisothermal calculations presented in this paper and in the 
isothermal calculations of Truelove et al. (1997), reflnement to Level 2 and 
above is controlled by application of the Jeans condition. This criterion can 
be augmented or replaced by others as necessary. The Jeans condition used 
to make the cell-by-cell refinement decision is a physically motivated reso- 
lution constraint discussed by Truelove et al. (1997). Jeans (Jeans 1902 and 
Jeans 1928) analyzed the linearized equations of 1-D isothermal GHD for 
a medium of infinite extent and found that perturbations on scales larger 



than the Jeans length, Aj = ( ^ ) , are unstable. Thermal pressure can- 



not resist the self-gravity of a perturbation larger than Aj, and runaway 
collapse results. Truelove et al. (1997) showed that the errors generated by 
numerical GHD solvers can act as unstable perturbations to the flow. In a 
simulation with variable resolution, cell-scale errors introduced in regions 
of coarser resolution can be advected to regions of flner resolution, aff^ording 
these errors the opportunity to grow. The unstable collapse of numerical 
perturbations can lead to substantial fragments, a process termed artificial 
fragmentation^ which can be avoided if Aj is sufficiently resolved. Defining 
the Jeans number resolution of Aj. Defining the Jeans number J = yj, 
Truelove et al. (1997) found keeping J < 0.25 avoided artificial fragmen- 
tation in isothermal evolution of a collapse spanning 7 decades of density, 
the approximate range separating typical molecular cloud cores from non- 
isothermal protostellar fragments. The constraint that \j be resolved is 
termed the Jeans condition. Although it has been shown to be necessary 
only for isothermal evolution, it is reasonable to expect that it is necessary 
(although not necessarily sufficient) for nonisothermal collapse as well. 

As a side effect of confining cell-sized perturbations to a length scale at 
which they are thermally damped, resolution of Aj also ensures that gradi- 
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ents developed in isothermal flow by gravity are well resolved. Formation of 
structure on scales of Aj and larger is a general feature of self-gravitational 
isothermal GHD flow since smaller fluctuations are damped but larger ones 
collapse. For example, in the self-similar solution for isothermal cylindrical 
collapse (Inutsuka & Miyama, 1992), the radial scale height of the cylinder 
is Aj. Lack of resolution of gradients within simulated flow triggers the in- 
jection of artificial viscosity, which is generally intended to be introduced 
in small amounts only for numerical stability (in our Godunov scheme) or 
shock mediation (in many common non-Godunov schemes). Introducing ex- 
cess amounts of artificial viscosity renders the problem solved different from 
the inviscid problem posed. Continuous resolution of Aj, however, keeps the 
flow inviscid and prevents artificial slowing of gravitational collapse. It is 
important to note, however, that the Jeans condition is a necessary but 
not, in general, sufficient condition to ensure convergence. 

The Jeans condition on Ax fundamentally differs from the Courant con- 
dition on At, although at first the two conditions might appear analogous. 
The Courant condition arises from a modal stability analysis of finite differ- 
ence equations derived from the Euler partial differential equations (PDEs) 
of hydrodynamics (see, e.g., Richtmyer &; Morton 1967). It is entirely a 
consequence of the finite differencing and has no physical counterpart in 
the PDEs. A modal stability analysis of finite difference equations derived 
from the GHD PDEs docs not yield the Jeans condition, but rather a gen- 
eralized Courant condition that includes effects of gravity (Truelove et al. 
1998). 

3. APPLICATIONS TO ROTATING MOLECULAR CLOUDS 

In this section we briefly describe recent high-resolution results wc have 
obtained with AMR for the gravitational collapse of nonisothermal, ro- 
tating, uniformly dense clouds (Fisher, Klein, & McKee 1998). Isothermal 
centrally condensed clouds (Truelove et al, 1997) and isothermal uniform 
clouds (Truelove et al, 1998) have been previously discussed and we refer 
the reader to these papers. 

The uniform cloud is an idealized model of an astrophysical cloud but 
is very useful as a first approximation. The initial isothermal configuration 
of the uniform cloud may be completely parameterized by the dimensionless 

energy ratios a = i^'thcrmal/ 1 ^gravitational I and/3f7 = £'rotational/|-£'gravitational| 

Burkert & Bodenehimer 1993 studied a uniform cloud, with M = IMq 
and R = 5 X 10^^ cm, which give po = 10~^^-^ g cm~^. Its energy ratios 
are a = 0.26 and f3ci = 0.16. The isothermal uniform cloud thus has Cg = 
0.167 km s~^ and a rotation rate Q = 7.2 x 10~^^ rad s~^. The cloud is 
perturbed hy p —>■ px [1-1-0.1 cos(20)], a 10% m = 2 mode seed perturbation. 
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This cloud begins with Aj = 1.17R. We have used an initial resolution of 
R32 and dynamically refined so as to ensure Jmax = 0.25. In the initial 
configuration, the fiducial radius beyond which cloud gravity dominates 
perturbation pressure is rcrit = 0.14i?, so that with R22 we still linearly 
resolve this scale with more than 4 cells. 

3.1. NON-ISOTHERMAL UNIFORM CLOUDS 

In recent work using high resolution AMR (Fisher, Klein, & McKee 1998), 
we followed the collapse of an initially rigidly rotating, uniform isothermal 
cloud. We used a two-component barotropic equation of state which makes 
the transition from isothermal to polytropic in a smooth fashion. The initial 
conditions are identical to the Burkert Sz Bodenehimer 1993 isothermal 
uniform cloud perturbed by a 10% m = 2 mode. In the following we discuss 
the evolution of this nonisothcrmal collapse. While a similar calculation 
has been recently studied by Bate&Burkert (1997), to our knowledge, none 
of the calculations in the literature to date have been able to follow the 
subsequent evolution of the fragments over dynamical (orbital) timescales, 
while still adhering to the Jeans criterion. However, to be fully certain of 
the reality of the fragmentation, it is still necessary to do a convergence 
study by decreasing the Jeans number (Truelove et al. 1997), which we 
have not yet done for this calculation. 

Figure 1 represents a 2-D slice of log density through the equatorial 
plane t = lAl x 10^^ s after start of the collapse. The cloud has initially 
collapsed to an isothermal disk, and an elongated filamentary bar with the 
first signs of fragmentation apparent. A strong isothermal shock above the 
plane of the disk is soon established. After t = 1.46 x 10^^ s (Figure 2) the 
isothermal bar becomes optically thick, and the accretion flow onto the bar 
is arrested, resulting in the growth of non-axisymmetric perturbations in the 
bar. Fragmentation in the bar results in the formation of binary spherical 
cores connected by a prominent bar. The core-bar system is embeddded in a 
two-armed spiral, derivative of the original m = 2 perturbation. The binary 
separation decreases as the cores increase their mass by direct accretion 
from the low-angular momentum material of the bar, eventually leading to 
the dissipation of the bar. Figure 3 shows a snapshot of log density about 
0.5 rotation periods later at t = 1.51 x 10^^ s. Protostellar disks have now 
formed around the cores with the cores at their closest orbital separation 
~ 44 AU. The disks accrete gas directly from the long spiral arms. Each 
core has O.O8M0 with a radius ~ 10 AU. The masses and radii of these 
first collapse cores are in good agreement with the detailed 1-D radiation- 
hydrodynamic calculations of Masunaga, Miyama, & Inutsuka (1998). The 
core accretion luminosity ~ 2 x 10^^ ergs" ^ obtained in the 3-D collapse 
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calculations is considerably less than that found in the 1-D simulations, 
~ 4 X lO^^ergs"^, due to the large angular momentum barrier in the 3-D 
simulations, resulting in slower accretion onto the protostellar cores. 

It is important to point out that the initial rapid growth of the cores 
is due to accretion of matter from the bar. As the bar dissipates, the cores 
have ended their first phase of growth and proceed to grow more slowly by 
direct accretion from the surrounding protostellar disks. At t = 1.52 x 10^^ 
s (Figure 4) the cores begin to separate and the surrounding disks become 
morphologically distinct. The protostellar disks and cores become a fully 
detached binary by t = 1.56 x 10^^ s (Figure 5) with the disks attached to 
the long spiral. The scale of the disks, of order 100 AU, is consistent with 
observations of gaseous disks surrounding single T-Tauri stars (Sargent & 
Beckwith 1994) and debris disks surrounding systems such as Pictoris 
(Artymowicz 1997). When we stopped our calculations, at t = 1.6 x 10-*^^ 
s(Figure 6) the binary protostellar disks/core system has moved back into 
the surrounding spiral arms. The cores have 20% of the mass of the initial 
cloud and the arms have 27% of mass of the initial cloud at this time. The 
protostellar disks comprise about 2% of the cloud mass. 

4. CONCLUSIONS 

In this paper we have described a powerful new method for numerical so- 
lution of 3-D self-gravitational hydrodynamic problems. This method com- 
bines a Godunov hydrodynamics integrator with a multigrid gravity solver 
in an adaptive mesh refinement framework. Guided by the Jeans condition 
for isothermal problems, and the use of Richardson extrapolation, AMR ef- 
ficiently expends computational resources only when and where the features 
of the flow demand it. We presented results of this methodology applied to 
the collapse and fragmentation of nonisothermal molecular clouds that are 
both initially uniform and centrally condensed, evolved over 8-9 decades 
of collapse and subject to initial m = 2 mode perturbations. We found bi- 
nary fragmentation results in good agreement with published calculations 
after evolution over the first half of this logarithmic range. By automati- 
cally inserting ever finer grids to maintain resolution of the Jeans length, 
these calcTilations did not generate the substantial artificial viscosity found 
in published calculations that used fixed finest resolutions. We have shown 
the first well-resolved formation sequence of protostellar disks surrounding 
a pair of binary cores, and discussed the role of the bar connecting the 
binary pair in the accretion onto the cores. The disks formed are consistent 
with observations of gaseous disks surrounding single T-Tauri stars and 
debris disks surrounding systems such as /3 Pictoris. 
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Uniform cloud with 10% m — 2 perturbation. 
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Global log p^^^ = -10.055221 I = 3.7159824 km s"' 

Resolution: 240 x 240 cells; Ax = 2.441 0000e+ 1 3 cm 
Plot radius: 0.058584100 R 
Run 333, Step 2300 

Cells 3976-4215 x 3976-4215 @ 4096 
chi_rho=100 
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Uniform cloud with 10% m — 2 perturbation. 
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XY slice of log density at equatorial plane 
Resolution level 4 of 4 
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log density 



Global log p^^^ = -9.9783973 I = 3.2023498 km s" 

Resolution: 240 x 240 cells; Ax = 2.441 0000e+ 1 3 cm 
Plot radius: 0.058584100 R 
Run 333, Step 3600 

Cells 3976-4215 x 3976-4215 @ 4096 

chi_rho=100 ^^^^H -I6.410000 
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Uniform cloud with 10% m — 2 perturbation. 
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XY slice of log density at equatorial plane 
Resolution level 3 of 4 
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log density 



Global log p^^^ = -9.9657727 I = 2.7501936 km s"' 

Resolution: 241 x 241 cells; Ax = 9.7659998e+ 1 3 cm 
Plot radius: 0.23537660 R 
Run 333, Step 4800 
Cells 904-1144 x 904-1144 @ 1024 
chi_rho=100 
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